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Abstract 

We present a new system for robustly performing Boolean operations on linear, 3D polyhedra. Our system is exact, 
meaning that all internal numeric predicates are exactly decided in the sense of exact geometric computation. Our 
BSP-tree based system is 16-28 x faster at performing iterative computations than CGAL’s Nef Polyhedra based 
system, the current best practice in robust Boolean operations, while being only twice as slow as the non-robust 
modeler Maya. Meanwhile, we achieve a much smaller substrate of geometric subroutines than previous work, 
comprised of only 4 predicates, a convex polygon constructor, and a convex polygon splitting routine. The use of 
a BSP-tree based Boolean algorithm atop this substrate allows us to explicitly handle all geometric degeneracies 
without treating a large number of cases. 

Categories and Subject Descriptors (according to ACM CCS): 1.3.5 [Computer Graphics]: Computing 
Methodologies—Computational Geometry and Object Modeling; 


1. Introduction 

Despite a long history of robustness issues. Boolean set op¬ 
erations (aka. constructive solid geometry or CSG for short) 
remain a popular choice in most 3d modeling systems avail¬ 
able. The operation conforms nicely to architectural concep¬ 
tions of positive and negative space in building massing, the 
machining of mechanical parts in CAD/CAM, and general 
purpose “sculpting” operations of adding and removing vir¬ 
tual stuff from a shape, to mention a few applications. We 
built the system described here to address such sculpting ap¬ 
plications, for which we found existing Boolean operations 
lacking in speed and robustness, particularly when perform¬ 
ing iterated operations, such as repeated gouging of an ob¬ 
ject with a tool. Sculpting is a particularly challenging appli¬ 
cation, since even moderate success requires fast(real-time) 
and completely robust(zero failure) solutions in order to be 
usable. We demonstrate that Boolean operations on polyhe¬ 
dra can be made fast and robust enough to be a viable com¬ 
ponent of such applications, contrary to the expectations of 
many in the computer graphics community. 

One of the first and biggest motivations for research into 
robust geometry was the not uncommon experience of see¬ 
ing a system crash in response to attempting to perform 
Boolean operations. However, despite 20 years of research 
on the problem, state of the art robust Boolean operations 
[HK05] may run up to 20 times slower than non-robust op- 
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erations. In our experiments(§4), these performance char¬ 
acterizations are actually optimistic. We see cases where 
the state of the art (CGAL) takes more than 50 times as 
long as non-robust Booleans. When we spoke to practition¬ 
ers who implement Boolean operations in entertainment and 
architectural design modelers, none were aware of any re¬ 
search into geometric robustness. They believed that there 
are no solutions to their robustness problems short of using 
Matlab/Mathematica-style big number computation, which 
they deem prohibitively expensive. Given that the aforemen¬ 
tioned Hachenberger et al. paper [HK05] constitutes the only 
experiment with more than a single test in the last 20 years, 
it is easy to see how this impression persists. 

Contribution We demonstrate a system capable of comput¬ 
ing completely robust (i.e. exact §2.1) linear Boolean opera¬ 
tions only twice as slowly as the non-robust (i.e. fragile §2.1) 
modeler Maya, providing the first empirical evidence that 
geometric robustness techniques are actually practical for 
computing Boolean operations. This constitutes the second 
and largest comparative test between non-robust commer¬ 
cial Booleans and robust Boolean systems ever conducted 
and published as far as we are aware. We accomplish this 
using a novel synthesis(§3) of techniques mostly 15 — 20 
years old(§2). In addition to a number of small but important 
changes to these techniques, we present a new plane-based 
convex polygon splitting algorithm(§3.2). 
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Figure 1: This bunny was sculpted out of a block using 
~ 9000 subtractions of dodecahedra. Our system computed 
these subtractions in approximately 3 hours at a rate of 
about 1 Hz- Using extrapolations from smaller experiments, 
this computation would have taken CGAL at least 3 days. 
Maya isn’t even able to perform this computation, crashing 
within the first couple hundred operations. 


We describe input/output methods(§3.4) for conversion 
between point-based polygonal meshes and plane-based bi¬ 
nary space partitioning trees, but I/O is not our focus. 
We make no strong guarantees about the accuracy of our 
conversions(e.g. topology preservation), but can perform 
the conversion arithmetic accurately to within machine ep¬ 
silon. Furthermore the current system performs robustly (no 
crashes/failures) for both plane-based and point-based in¬ 
puts, additionally ensuring 100% accurate results for inputs 
already represented by planes, such as when sculpting with 
primitive polyhedra—our original motivation. 

2. Background 
2.1. Robustness 

Following Yap [Yap04], “nonrobustness refers to qualita¬ 
tive or catastrophic failures in geometric algorithms aris¬ 
ing from numerical errors.” That is, geometric robustness 
is not the same as precise numerics. Numeric precision for 
its own sake is usually irrelevant, so long as the results 
aren’t grossly inaccurate. However, small numeric errors 
will sometimes lead to catastrophic program failures (incon¬ 
sistencies) or grossly incorrect results. Geometric robustness 
seeks to avoid these byproducts of numeric errors without in¬ 
curring the cost of fully precise numerics. Thus, geometric 


robustness is and has always been about simultaneous speed 
and robustness. Either goal is easy to achieve in the absence 
of the other. 

Numeric computations in a geometric program may be 
classified as either predicates or constructions [She06], 
Predicates make two-way or three-way decisions based on 
numeric coordinates of geometric objects, while construc¬ 
tions compute the coordinates for new geometric objects 
from the known coordinates of existing geometric objects. 
In general, algorithms that restrict themselves only to pred¬ 
icates are easier to make robust and fast than algorithms 
which use constructions. By restricting to predicates, each 
of which makes a decision based on a few numbers, the 
depth of any arithmetic expression has an a priori constant 
bound. This bound allows for static filters [She97, FVW96] 
that produce fast predicates that are exact, meaning they al¬ 
ways give the correct branching answer. By way of contrast, 
algorithms that use constructions allow for the construction 
of arbitrary depth arithmetical expressions, incurring sig¬ 
nificant speed penalties to guarantee exact answers. In this 
work, we make use of Sugihara and Iri’s observation [SI89] 
that linear Booleans require no constructions if representa¬ 
tions are based on planes and not points. 

Our solution is exact in that it is unconditionally robust 
given consistent input. It is also a fixed-precision tech¬ 
nique since we can, by the use of plane-based representa¬ 
tions, avoid constructions and thus provide an a priori con¬ 
stant bound on the depth of any arithmetic trees used and 
thus on the representational complexity required. This is 
in contrast to arbitrary-precision techniques which have 
no such bound, generally because they use constructions 
[BMP93, GHH*03]. As we will demonstrate (§4) fixed- 
precision techniques are much faster than arbitrary-precision 
techniques. Devillers and Pion [DP03] provide corroborat¬ 
ing evidence in their Delaunay triangulation experiments. 

There are a variety of non-exact approaches that provide 
weaker notions of robustness. Well-known examples include 
epsilon-tweaking [NAT90, LTH86], in which heuristic tol¬ 
erance parameters are adjusted in hope of successful exe¬ 
cution. Industrial systems for CAD/CAM and solid model¬ 
ing generally employ this technique along with undo oper¬ 
ations to allow the user to work around operations that fail, 
but it is a fragile system in that it provides no guarantee 
against failure of an operation. Interval geometry techniques 
[Seg90, Bru91], which represent geometric primitives with 
local tolerances within guaranteed intervals, are more princi¬ 
pled. They are quasi-robust in that they provide guaranteed 
bounds on the imprecision of the representation, but these 
intervals can become arbitrarily wide, eventually collapsing 
all the geometry to a point within a huge interval. As with 
epsilon-tweaking, systems based on this approach can be 
quite useful if users are willing to work around these prob¬ 
lems. However an efficient exact system is preferable (and 
sometimes necessary) in order to avoid user workarounds. 

©2009 The Authors) 
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General Position and Geometric Degeneracies 

In addition to the aforementioned numeric issues, many al¬ 
gorithms rely on “general position” assumptions, such as, 
“In general, two surfaces (2d) intersect in some number of 
curves (Id).”. Local arrangements of geometry that violate 
these assumptions are known as (geometric) degeneracies, or 
degenerate arrangements. One approach to geometric degen¬ 
eracies is to not make general position assumptions—we do 
this. However, for some algorithms this approach leads to a 
large, unwieldy explosion of degenerate cases. Thus geome¬ 
try is often perturbed into general position. Done naively, (by 
directly modifying coordinates) perturbations produce fur¬ 
ther complications. Therefore, robust perturbation [Sei94] is 
symbolic, further necessitating the use of exactness. 

2.2. Booleans 

We (along with Maya and industry) follow Requicha’s regu¬ 
larized sets formalism [Req77]. CGAL follows the Nef for¬ 
malization [Nef78]. 

There are many specific approaches to handling Booleans 
in the literature, but since we have to interoperate with sur¬ 
face modelers, we are concerned with those which support 
boundary evaluation and thus can output meshes. B-rep al¬ 
gorithms [HHK89, LTH86] are the most common of these, 
and they share an algorithmic template with the more recent 
cell complex algorithms [GHH*03, HK05] . This template 
is: 1) If A and B are the boundaries of two objects whose 
union, difference or intersection we would like to compute, 
find the intersection of A and B, thus dividing each surface 
into two components, one inside and one outside the other 
surface. 2) Select the appropriate component of each sur¬ 
face, and 3) stitch these together to form the correct out¬ 
put. This apparent simplicity belies the large number of spe¬ 
cial cases that result from the various ways the two objects 
can align [Hof89]. Resulting systems can be quite complex, 
making it difficult to ascertain whether all cases have been 
handled. This case analysis can be reduced by the use of per¬ 
turbation techniques(§2.1) [For97], but we have yet to come 
across any publication that explicitly presents a complete ge¬ 
ometric substrate for a B-rep algorithm. 

CSG trees encode objects as a Boolean combination 
of primitive objects (e.g. (A U (B n C)) — D). Comput¬ 
ing Boolean operations on CSG trees is often straightfor¬ 
ward, the real problem is producing the resulting bound¬ 
ary [RV85], This can be accomplished by reduction to B- 
rep algorithms or by algorithms that operate globally on 
the tree and exploit properties of a restricted set of primi¬ 
tives [Bru91, BR96]. However, these methods are inappro¬ 
priate for iterative, and thus interactive applications. 

BSP trees afford an alternative to B-rep algorithms that 
avoid their concomitant case explosion by explicitly han¬ 
dling all degenerate configurations of geometry [TN87, 
NAT90]. These have been demonstrated to have suitable 
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performance for interactive volumetric sculpting [Nay90], 
and recent work [LDS09] suggests further possible per¬ 
formance improvements. Unfortunately, all of these ap¬ 
proaches are fragile and therefore subject to failures and user 
workarounds. We will detail this approach shortly(§3.3). 

Aside from the BSP tree work, the closest precursor to 
our system is that of Fortune [For97]. Like us, Fortune em¬ 
ployed both Sugihara and Iri’s observation on planes and 
fixed-precision techniques for robustness. Beyond this, we 
have very little in common. We use different predicates, a 
completely different representation: BSP trees instead of B- 
reps, and we do not rely on symbolic perturbations, instead 
choosing to handle all degenerate arrangements explicitly 
(and with a simpler substrate). These differences are all the 
result of key design decisions toward our goal of combined 
robustness and performance. 

Sampled volumes offer another alternative—gridded 
discretization—often advocated for sculpting contexts 
[GH91], This approach was refined by Ferley, Cani and 
Gascuel [FCG00], and reached its zenith with Perry and 
Frisken’s Kizamu [PF01], capable of handling up to 12-ply 
deep octrees. Two major shortcomings of Kizamu dissuaded 
us from pursuing volumetric approaches further: (a) render¬ 
ing becomes prohibitively expensive as 12-play octrees are 
approached and (b) discretization errors, more deleterious to 
perceived visual quality than roundoff, are still present. 

3. System Design 

Although we began by trying to build a Boolean sculpting 
system based on Thibault and Naylor’s BSP tree algorithms 
[TN87,NAT90], it is, perhaps, easier to motivate and explain 
our design in reverse; that is, bottom up. As Sugihara and 
Iri pointed out, one must use planes, not points as the funda¬ 
mental primitive in order to make use of fixed precision tech¬ 
niques (§2.1). We proceed from this decision upwards. Nat¬ 
urally, our numeric predicates must all be phrased as predi¬ 
cates about planes. Likewise, the geometric definitions and 
subroutines must now bootstrap off of these predicates and 
thus be completely phrased in terms of plane manipulations. 
Finally, binary space partitions are employed because they 
are (a) plane-based, and (b) allow for economy in the geo¬ 
metric substrate. 

3.1. Numeric Substrate 

A plane, in our system, is taken to be a quadruple of float¬ 
ing point numbers, (a,b,c,d) interpreted as the coefficients 
of a plane equation. On occasion we will have cause to re¬ 
fer to “points”, by which we do not mean coordinate triples. 
Rather, we take a point to be a triple of planes ( p,q,r ) whose 
mutual intersection is the point in question. We do not use 
or refer to lines. We provide exactly four arithmetical pred¬ 
icates concerning planes and their relative positioning in 
space [YN97,Sto87], 
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coincidence. Two planes p,q are coincident, if and only if 
the determinants of all 2 x 2 minors of the following matrix 


T Pa Pb Pc Pd ] 

[ qa qb q c qd \ 


Construction of a Convex Polygon from a Plane Given a 
plane h, this operation constructs a convex polygon repre¬ 
senting h clipped by a “very large box.” In this way, the out¬ 
put polygon serves as a stand in for the infinite extent plane. 
A consistent axis aligned “very large box” is used for all 
calls to the constructor and consists of the volume of space 
bounded by the planes x + , x ~, y + , y~, z + , z ~. 


coincident orientation. Furthermore if p and q are known 
to be coincident, then p is similarly oriented to q if and only 
if all products p a q a , Pcqc, Pdqd are non-negative. 

point validity. A point (p, q, r) is well-defined (i.e., valid) 
if and only if the following determinant is non-zero: 


I Pa Pb Pc I 

qa qb qc 

| r a r b r c \ 

orientation. Given that the point ( p,q,r) is valid, it lies 
behind, on, or in-front of the plane s if and only if the fol¬ 
lowing expression is negative, zero, or positive, respectively: 


I Pa Pb Pc 

qa qb qc 

| r a r b r c 


Pa Pb Pc Pd 

qa qb qc qd 

r a r b r c r d 

Sa s b Sc S d 


These predicates are implemented as static filtered 
floating-point predicates in the style of Shewchuk [She97]. 
We make a few deviations from Shewchuk’s pattern. First, 
our predicates only take single-precision floating point input, 
but perform computation using double precision. This allows 
us to guarantee the absence of exponent overflow/underflow, 
and allows us to tighten our static error bounds. Second, 
rather than using a cofactor expansion, we compute the 4 x 4 
determinant as a dot product between two lines in Pliicker 
coordinate form. This makes the arithmetic tree shallower, 
resulting again in tighter static error bounds. Finally, we only 
use single stage filters, for simplicity’s sake. 


3.2. Geometric Substrate 

In the preceding section on the numeric substrate we de¬ 
fined planes as primitives, points as triples of planes and 4 
predicates operating on these planes. We will now define a 
convex polygon type along with a constructor and splitting 
routine. These 2 subroutines constitute the entire geometric 
substrate. They are used to support sub-hyperplane construc¬ 
tion, tree splitting, and input/output conversions(§3.3, §3.4). 

convex polygon. We take a convex polygon (In this paper, 
‘polygon’ always stands for ‘convex polygon’) to be a plane 
of support v along with a list of bounding planes {h;}; e g n 
listed in counter clockwise order. The vertices of this poly¬ 
gon are then given by v,- = (s, b;_ i, b\). 


In order to form the resulting polygon the dominant com¬ 
ponent of h’s normal is found. Without loss of generality, let 
it be h z . Then a polygon P is formed combinatorially with 
h as support and x + ,x~ ,y + ,y~ as bounds. The order of the 
bounds may be determined by inspecting the sign of h z . Fi¬ 
nally, P is clipped by z + and z _ using the (following) poly¬ 
gon splitting routine to ensure that the entire polygon lies 
within the “very large box.” 

Splitting a Convex Polygon by a Plane Since polygon 
splitting may be viewed as two successive and complemen¬ 
tary instances of polygon clipping, we will present a polygon 
clipping algorithm for simplicity. Given a convex polygon 
(.?. {A' jigzJ, and clipping plane h, we wish to output the 
part of the polygon lying to the positive side of h, or indicate 
that there is no such part. 

Our algorithm is similar to Sutherland-Hodgman style 
polygon clippers, where one decides at each step “do or do 
not output the current vertex”, except we decide “do or do 
not output the current bounding plane (edge).” Analogously, 
rather than deciding if and when to insert a crossing point 
into the output stream, we decide if and when to insert the 
splitting plane into the output stream. These changes are key 
to eliminating the construction of line-plane intersections, il¬ 
lustrating the “use planes, not points” principle. 

If i is coincident with h 

If s is similarly oriented to h 
Return (s, {b,} ieZ J 

Else 

Return nothing 

Else 

For i £ Z„ (in order), 

Output planes as specified by table lookup using: 

o(s,bi_2,bi-i,h) 

o(s,bi~ i ,bi,h) 

o(s,bi,b i+ l,h) 

Return (s , output) 

Algorithm: Clip (v, {b;};exj by h 

The algorithm proceeds by iterating through the poly¬ 
gon’s bounding planes and deciding whether to output each 
bounding plane £>,. In order to make this decision the two 
endpoints of h;, v,- = (j,b,_i ,bj) and v* + i = (s,bj,bj + 1) are 
tested against the clipping plane h. Additionally, the imme¬ 
diately preceding vertex, v,_ [ = (s,b,-_2,b;_i) is tested. As 
we proceed, the reason for testing this additional third vertex 
will manifest itself. 
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Figure 2: (left) coding example (center & right) These two 
cases are indistinguishable using only two vertices. 


Given that each vertex may lie behind (—), on (0), or 
in-front (+) of the clipping plane h, each triple of vertices 
v,-_i, v,-,v !+ i may lie in one of 27 possible arrangements rel¬ 
ative to the clipping plane h. Asa shorthand we use the fol¬ 
lowing coding scheme. Suppose that v,_i lies on (0) h, that 
v, lies behind (—) h, and that Vj+i lies in-front (+) of h. Then 
encode this arrangement with the sequence 0— h (Figure 2) 

For each of these 27 codes, we must decide whether to 
output the current bounding plane b; (signal B) or not (no 
signal); we must decide whether to output the clip plane h 
(signal FI) or not (no signal); and if both are to be output, 
then we must decide on the relative order of output (BFf or 
HB). To complete the procedure, the output planes (if there 
are any) are taken in order as the hounding plane list for the 
resulting clipped polygon. 

In deriving the correct output, one should be careful of the 
following two pitfalls. One must make sure that the H signal 
is given no more than once, and one must not output both the 
bounding plane b ,■ and the clip plane h in the case that bf s 
corresponding edge (v, to v* + i) lies in the clip plane h. We 
address the former pitfall by choosing to output h only when 
re-entering the positive halfspace. We address the latter pit- 
fall by omitting the output of the hounding plane in question 
bi, and by outputting the clip plane h as the polygon turns 
back into the positive halfspace. 


Table 1: Polygon Clipping Encodings (* is a wildcard) 


Input 

Output 

I 

0 

I 

o 

*++ 

B 

+ 0 + 

B 

*-+ 

HB 



00 + 

HB 





-0 + 

HB 



*+0 

B 

*00 

0 

*-0 

0 

*+- 

B 

*0- 

0 

*— 

0 


Upon close inspection of this table, you can see that if 
the +0-1- case output HB instead of B, then the entire table 
would reduce under symmetry to three or four cases, and 
appear even more similar to Sutherland and Hodgman style 
clippers. To understand why this asymmetry occurs, con¬ 
sider Figure 2, with two instances, one of a hexagon clipped 
along its diagonal and one of a triangle clipped at a single 
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vertex. Using only two vertices, we do not have sufficient 
local information to distinguish whether or not the clipping 
plane should be output. 


3.3. BSP Tree Set Operations 

Algorithms for Boolean set operations on polyhedra using 
BSP trees were presented by Thibault, Amanatides and Nay¬ 
lor [TN87,NAT90], This approach has the advantage over B- 
rep based approaches of handling degeneracies in a simpler 
and more uniform way, but it has been prone to numeric ro¬ 
bustness issues. We use a modification of this approach that 
is completely plane-based so as to allow for use of our effi¬ 
cient, robust substrate. We present a high-level overview of 
the BSP tree approach, followed by a description of our de¬ 
viations from the algorithms that are described in full detail 
in Thibault and Naylor [TN87] and Naylor et al. [NAT90], 



■ in - a + 

■° ut x° ut 

/\ 0UT 

/ d \ /\ 

IN OUT IN OUT 


Figure 3: BSP tree representing the blue object 



sub-hyperplane 


Region 


Boundary in this 
sub-hyperplane 


Figure 4: parts of a BSP tree internal node (cell diagram) 


How does a BSP tree represent an object? Using planes, 
a BSP tree recursively partitions space into (convex, possi¬ 
bly unbounded) cells. The leaves of the BSP tree correspond 
to indivisible cells in the partition. By simply coloring each 
of the leaf cells with one of the two labels, IN or OUT, we 
get a representation for the polyhedral set consisting of all 
of the IN colored cells. (Figure 3) Internal nodes of the BSP 
tree represent divisible cells of the partition. Each such cell 
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has an associated “region”, which is a set of halfspaces de¬ 
scended on a path to this node from the root. The intersection 
of these halfspaces defines this region, ie. the extent of the 
cell, geometrically. At the node of the tree, the partitioning 
(hyper)plane is stored, which along with the region allows 
for finding the sub-hyperplane, that part of the hyperplane 
lying in the node’s region. Additionally, positive and nega¬ 
tive child pointers are maintained. (Figure 4) 

How does a BSP tree represent the boundary of an ob¬ 
ject? The boundary is stored at the internal nodes of the BSP 
tree as an augmentation of the tree. Each internal node of the 
BSP tree stores the portion of the boundary contained within 
its associated sub-hyperplane, represented as a set of con¬ 
vex polygons. For us, these are plane-based polygons. For 
Thibault and Naylor, they were point-based polygons. 



How do you compute set operations, say A U B, between 
BSP trees? The two operand (binary space) partitions are 
“merged” into a new (binary space) partition (Fig 5). This 
new partition is found by repeatedly splitting one tree by the 
other’s root partition, thus reducing the set operation prob¬ 
lem to two smaller instances at each step. The leaf cells of 
this final partition are either IN or OUT in each operand ob¬ 
ject A and B. Thus, the resulting leaf is colored IN if it is IN 
A or B and colored OUT if it is OUT of both A and B (and 
similarly with different logic for intersection and difference). 

How are degeneracies dealt with? (i.e. How is the tree 
split?) Given a subtree T lying in a region (Fig 4), and a 
plane to split it with S, we split T’s sub-hyperplane ( T.shp ), 
and then recursively split T’s subtrees. In order to do so, we 
need only determine how T.shp and S IT T.region lie with 
respect to each other. Figure 6 shows the 7 possible arrange¬ 
ments that arise. Determining which of these 7 hold is a sim¬ 
ple matter of splitting T.shp by S and S fl T.region by T.hp. 

How do you compute the boundary of this new, resultant 
object AUB? The new boundary must be composed of some 
subset of the two old boundaries. Therefore, we can just clas¬ 
sify the old boundaries, using the new tree, to see which 
faces we should keep around. 



Figure 6: (a) T.shp and SC\T.region are coincident and 
either similarly or oppositely oriented, (b) T.shp and S IT 
T. region do not intersect, and are relatively oriented in one 
of four possible combinations, (c) T.shp and S DT.region 
intersect. 


We can implement most of these operations recursively, 
over the BSP tree and its subtrees. Furthermore, we can 
be very efficient by interleaving the computation of the 
new boundary and BSP tree together. By using information 
gained in the midst of the computation, we can avoid un¬ 
necessary work, e.g. the recursive algorithm can quickly de¬ 
termine that a whole subtree of A lies disjoint from B, and 
therefore the boundary stored in this subtree is not modified 
by the set operation. In this way, BSP trees naturally exploit 
spatial locality by equating it with the structural locality of 
the tree structure. 

Our particular implementation is a hybrid of techniques 
from Thibault and Naylor’s two papers [TN87,NAT90]. We 
take the “Merge” algorithm from the ’90 paper as our scaf¬ 
fold. Building onto that, we store the boundary as convex 
polygons at the internal nodes and maintain their consistency 
throughout the set operation, as advised. In a departure from 
Thibault and Naylor’s work, we do not worry about union¬ 
ing or gluing together boundary fragments. Instead, we allow 
for boundary polygons within a sub-hyperplane to overlap. 
If a non-redundant boundary is needed (say for output) then 
we reconstruct a boundary, directly from the tree, for output. 
This change gives us one less thing to worry about. From 
the ’87 paper, we take the BSP tree reduction techniques. 
We found tree reduction to be very important for saving both 
time and space. 


3.4. Interoperability: Input and Output 


Plane-based Representations 


Point- /1L Plane- 
based based 
Polygon r-^\Polygon 
Soup L n/ Soup 


BSP Trees 
(plane- 
based) 


The system we have described is designed to work within 
a nice sanitized world of objects represented as plane-based 
BSP trees. Unfortunately, no other operations or systems are 
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designed to work within this world. In order to accommodate 
these circumstances, we provide the following translations: 
between point and plane-based polygon soup, and then be¬ 
tween plane-based polygon soup and BSP trees. In this way, 
input and output to and from the plane-based BSP tree rep¬ 
resentation may be decomposed into 4 simpler transforma- 

point soup —> plane soup. Each point-based polygon is 
approximated by a plane-based polygon. First, a plane is fit 
to the vertices of the point-based polygon, to serve as the 
plane of support s of the new plane-based polygon. Using 
the plane to polygon construction (§3.2), an initial polygon 
is created. Then for each edge of the point-based polygon, a 
new bounding plane is fabricated, passing through that edge 
and orthogonal to the plane of support s. The initial polygon 
is then successively clipped by these bounding planes, using 
the polygon splitting routine (§3.2). 

plane soup —► BSP tree. Because we are not guaranteed 
that the resulting plane-based polygon soup is watertight, 
we assume the worst and implement a mesh repair algo¬ 
rithm. Murali and Funkhouser have devised a clever mesh 
repair scheme, based on BSP trees [MF97]. The mesh repair 
scheme builds a consistent BSP tree as an intermediary step, 
and so is perfectly suited to this role. Alternatively, if we 
know that the plane-based polygon soup is water-tight (for 
instance, if it was a boundary extracted from a BSP tree), 
then we can use a standard BSP tree building algorithm, such 
as the one described by Thibault and Naylor [TN87], 

BSP tree —* plane soup. Since the boundary is already 
stored in the BSP tree, all we need to do is traverse the tree 
and collect it. As a minor note, we may want to recompute 
the boundary of the tree once before doing so for reasons 
previously mentioned(§3.3). 

plane soup —> point soup. Since the vertices of a plane- 
based polygon are implicitly defined as the intersection of 3 
given planes (§3.2), we may convert to point-based output by 
computing/approximating the coordinates of these vertices. 

connectivity. Although we never use or record connec¬ 
tivity information (as one would when working with B-reps) 
many other programs expect to find connectivity information 
encoded in their input polygon soup. All of this information 
may be recovered by knowing precisely when two vertices of 
two different polygons are actually the same point in space. 
Since any plane based polygon soup output from a BSP tree 
is guaranteed to be valid, we can recover this vertex inci¬ 
dence information from the plane soup using the orientation 
predicate. Given two valid points ( a,b,c ) and (x,y,z) (de¬ 
fined as triples of planes), ( a,b,c ) is incident to (x y, z) if 
and only if (a, b, c) lies on all three planes x, y, and z. 

accuracy. For our purposes (sculpting) accuracy was only 
important (a) up to visual inspection and (b) in order to 
achieve robustness. Towards this end we make a distinc¬ 
tion between the accuracy of conversion (the consistency 
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between the input and output of a given conversion step) 
and validity of output (the self-consistency of a represen¬ 
tation). Note that the proposed I/O path always ensures va¬ 
lidity. Therefore the program will not crash or fail. However, 
the I/O steps do not always guarantee 100% accuracy. In par¬ 
ticular, conversion to and from point and plane based poly¬ 
gon soup is computed inexactly, making it subject to round¬ 
off error. Numeric accuracy is only lost in these conversions. 
(This loss may be mitigated to machine epsilon in both input 
and output sides using Priest’s techniques [Pri91], on which 
Shewchuk’s predicates are based.) Additionally, the pertur¬ 
bation caused by roundoff during input will break every non- 
trihedral vertex (with degree k > 3) into k — 2 nearly coin¬ 
cident, but distinct vertices, potentially introducing microp¬ 
olygons and/or splinter-like polygons. Because we are only 
concerned with accuracy up to visual inspection, we can run 
a mesh simplification program on our output to collapse split 
vertices and eliminate degenerate polygons. Since this mesh 
simplification only relies on valid connectivity data to oper¬ 
ate robustly, no new robustness issues are introduced. 

4. Experiments 

4.1. Setup 

We ran four tests on our system (BSP), CGAL’s [CGA08] 
Nef polyhedra, and Maya’s [Aut09] Booleans. Under our 
taxonomy (§2) BSP uses an exact (fixed precision) BSP tree 
approach, CGAL an exact (arbitrary-precision) B-rep (tech¬ 
nically Nef polyhedron) approach, and Maya a fragile B-rep 
approach. We ran all tests on a 1.83 Intel Core Duo (32bit) 
Macbook Pro. We compiled CGAL with static libraries, 
linked using gcc and options -02 -DNDEBUG. CGAL ker¬ 
nels were Exact_predicates_exact_constructions. We used 
A trial version of Maya 2009. No paging, etc. occurred. 

We report times for performing Boolean operations, 
whose definition varies from system to system. For BSP, this 
consists of computing the resulting tree with its boundary 
embedded(§3.3). For CGAL, this consists of computing the 
resulting Nef polyhedron. However, some extra computation 
(beyond lookup) is then necessary to compute the boundary. 
We did not measure this time since we felt it would unduly 
punish the design decision to use Nef Polyhedra. For Maya, 
MEL scripting was used to measure the time to execute a 
Boolean command. Input conversion times (only applicable 
to octoball) are provided separately from Boolean times. 

4.2. Tests and Results 

This battery of tests was chosen to give a good coverage of 
potential Boolean use in a sculpting system from single op¬ 
erations between two sizable meshes (Octoball) to iterated 
applications of small tool shapes (Random Boxes, Sculpt 
Bunny). We also included the heatsink test to check for worst 
case performance when the output is 0(n 2 ) in the size of the 


G. Bernstein & D. Fussell / Fast, Exact, Linear Booleans 


Octoball Test The Octoball test computes a single set oper¬ 
ation: a lumpy ball volume minus an octopus volume (Figure 
7). The resulting ball is left with many tunnels and a few in¬ 
ternal bubbles (octopus eyes). In addition, this test was use¬ 
ful as an “organic” rather than synthetic test, since the mod¬ 
els were found rather than fashioned for this purpose. We 
scaled this test by applying mesh simplification to the octo¬ 
pus and ball models, yielding 10 copies with 0.1,0.2,..., 1.0 
times as many faces as the original meshes, which had 11444 
faces (ball) and 33872 faces (octopus) respectively. 

No systems failed on any instance of this test. BSP ran 
roughly 1.5x slower than Maya and about 3-6x faster than 
CGAL. The input routines for BSP resulted in trees contain¬ 
ing 4x as many nodes as faces in the input model, due to mi¬ 
croscopic and near degenerate polygons formed by point to 
plane conversion and the subsequent mesh repair algorithm. 
However, BSP still maintained a memory footprint smaller 
than CGAL: BSP took 30MB for the 0.1 x face meshes com¬ 
pared to CGAL’s 50MB footprint. For the l.Ox face meshes 
BSP took 250MB vs. CGAL’s 270MB. (Memory usage was 
measured using process monitoring.) Furthermore, BSP took 
20 seconds and 6 minutes respectively for the 0.1 x and l.Ox 
conversions as compared to CGAL’s 6 minutes and 2 hours. 


Figure 7: (left) Octoball: the octopus is subtracted from the 
bumpy ball in the test. They are displayed unioned for visual 
reference; (right) Heatsink: the test result is displayed for 
n = 15; 225 rods 


Heatsink Test The heatsink test consists of computing a 
single intersection operation between n slabs arrayed along 
the x axis (vertical slabs) and n slabs arrayed along the y axis 
(horizontal slabs), resulting in n 2 rods (Figure 7). This test 
induces the worst case 0(n 2 ) behavior. 

Maya produced incorrect results on this test. BSP was 
much faster than either Maya or CGAL (asymptotically so 
and by an average constant factor of 50 x over CGAL). We 
do not know why Maya’s performance on heatsink appears 
to be exponential (see table), and we are loath to conjecture 
without further information. 

Random Box Test The Random Boxes test consists of n it¬ 
erative, varying unions and subtractions of randomly chosen 
rectangular prisms. The scene begins with the 0th box and 


each ith box thereafter is unioned into the scene unless i is 
divisible by 5 or 7, in which case it is subtracted. Box width 
is bounded between 1.0 and 3.0 in each dimension, positions 
are chosen randomly throughout a 10.0 X 10.0 x 10.0 cubic 
volume, after which a random rotation is applied. This test 
provides a synthetic stress test for iterated Boolean opera¬ 
tions, as one would use in sculpting. The same set of random 
boxes was used across all three tested systems. 

Maya failed this test gracefully after about 90 operations. 
By periodically rebuilding the tree (from plane-based poly¬ 
gon soup) every 20 operations, we were able to get a 2-3 x 
speedup on BSP. Furthermore, this dramatically reduced tree 
size, maintaing it beneath 20,000 nodes, rather than it pro¬ 
ceeding upwards of 100,000. As a result, BSP was much 
faster than CGAL (22-28 x). 

Sculpt Bunny Test The Sculpt Bunny test consists of au¬ 
tomatically sculpting the Stanford bunny out of a block of 
“clay” by repeated application of boolean difference oper¬ 
ations to subtract out smaller chunks. The chunks to be re¬ 
moved (dodecahedra) are placed by sampling the bunny vol¬ 
ume with a regular grid and choosing a subset of the samples 
outside the bunny. Radii (dodecahedra scaling) are associ¬ 
ated with each of these sample points so that the subtracted 
volume covers all of the samples outside of the bunny. This 
procedure was used to generate our Stanford Bunny model 
(Figure 1), using an originally 80 3 sample grid. For this test 
we used a 10 3 sample grid, ultimately resulting in 250 do¬ 
decahedra to subtract. This same set of dodecahedra was 
used across all three systems. 

Maya successfully completed these 250 operations. 
Therefore we tried a larger 20 3 resolution yielding about 
850 operations. On this test set, Maya gracefully failed after 
25 operations. Again, we used tree rebuilding to accelerate 
BSP, yielding a 2x speedup. With tree rebuilding, BSP was 
only 2x slower than Maya. BSP was again much faster than 
CGAL (16-17X). 

4.3. Analysis 

As a general trend, BSP was much faster (16-28 x iterated 
tests, 3-6x octoball) than CGAL, our exact modeler, yet 
only 1.5-2x slower than Maya, our fragile modeler. Given 
that our code has not been optimized, we believe this pro¬ 
vides strong evidence that we can build exact Booleans that 
run as fast as fragile Booleans. We have also definitively 
shown the excessive cost of arbitrary-precision exactness 
techniques for Booleans, even in the case of a single Boolean 
operation (octoball, heatsink). 

5. Conclusion 
5.1. Discussion 

We presented our system design(§3) as if we selected the 
constituent parts independently: one algorithmic approach, 
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one substrate and one robustness approach. However, these 
choices are highly dependent on one another. In order to 
get robustness on the cheap, we use fixed-precision robust¬ 
ness, not arbitrary-precision. To do so we must avoid con¬ 
structions, which requires using planes not points in our 
geometric substrate. Thus the interface provided to the al¬ 
gorithm must be plane-based, making BSP tree algorithms 
ideal since BSP trees are completely based on planes that 
partition space. In return the BSP tree algorithm requires an 
unusually concise substrate to handle all degeneracies. We 
collect these disparate reasons under the unifying maxim, 
“use planes, not points,” a fundamental insight into the na¬ 
ture of computing Booleans. 

Furthermore, swapping out a single component is not gen¬ 
erally possible. Without exact arithmetic, BSP trees become 
intolerably fragile. Even plane-based B-rep algorithms yield 
complicated substrates. Demanding points forces the use of 
arbitrary-precision arithmetic. 

5.2. Applications and Limitations 

Although our system was designed for sculpting, it seems 
reasonable to use in more traditional 3d modelers to re¬ 
place fragile Booleans. However, such systems are likely to 
use B-reps, not BSP trees. The described input and output 
paths(§3.4) provide a rudimentary sketch of the representa¬ 
tion conversions needed, but there are two major limitations: 
the conversion is slow (a few seconds to a few minutes), and 
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the result may be verbose and ill-conditioned (microscopic 
and near degenerate polygons). Given these limitations, our 
system could be most profitably used in a distinct “BSP tree” 
mode, where users explicitly convert models into BSP trees, 
perform set operations, and then convert back, allowing am¬ 
ple time for conversions. 

Other limitations include restriction to linear geome¬ 
try and not handling rotations. Rotation arithmetic induces 
roundoffs, inadvertently altering the arrangement of geome¬ 
try. As a quick solution, it’s possible to convert a BSP tree 
to plane-based polygons, rotate, and then convert back into 
a BSP tree (using mesh repair). 

5.3. Future Work 

As demonstrated (Figure 1) our system can compute 
Boolean operations at approximately 1Hz for reasonably 
complicated objects. One avenue for future work is trying to 
improve performance until robust Booleans are truly inter¬ 
active for complex objects. We see further performance im¬ 
provements coming in two ways: managing model complex¬ 
ity and special casing operations. Iterated operations tend to 
increase model complexity, slowing down both future mod¬ 
eling operations and rendering. This complexity manifests 
in tree size (number of nodes) and number of polygons. 
What sort of (possibly lossy) tree simplification can we per¬ 
form while still supporting Boolean operations? How can we 
make operations faster? Can we exploit simplifying assump- 
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tions, such as: The tool object is small relative to the sculpted 
object, it is convex, or repeated applications are likely to hap¬ 
pen close together? 

Another avenue for future work is improving conversion 
between BSP trees and polygonal meshes. Again, making 
the problem more specific will help. Given a consistent 
point-based mesh, can we skip mesh repair on input and 
get a consistent BSP tree? Can we use partial conversions 
to safely and quickly support other operations like rotations 
or even free-form deformations? 
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